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ABSTRACT 


This thesis deals with the problem of Digital 
Interpolation, or the process of increasing the sampling 
rate of a sampled data signal. It is mainly concerned with 
investigating the possibility of improving the performance 
of interpolators for bandlimited signals regarding which 
no additional information is available. The approach 
suggested is to estimate the MEM spectr\jm and design optimum 
interpolators based on this estimate. Using this approach, 
extensive studies on a representative set of signals were 
carried out. A comparative evalu0,tion of the performance 
of the various interpolation schemes was made on the basis 
of minimm mean square error criteria. The degree of 
improvement resulting from the suggested method is shown 
to depeiid on the natiore of the spectrum. The thesis also 
includes a review of work in Digital Interpolation, 
various interpretations and approaches to the problem 
and useful computer programs used in the study. 



CHAPTER 1 


DTTRODUCTIOH 


In. digital signal processing there ib often a need 
to alter the sampling rate of signals. Digital Interpolation 

refers to the process of increasing the sampling frequencT- hy 
filling in more samples in between those available in the 
input data. A related process, one of decreasing the sampling 
rate, is referred to as decimation. In most applications, the 
tv>ro operations are employed in a complementary manner. 

At the outset, the distinction between interpolation 
as is to be understood in digital processing and its usual 
concept in numerical analysis should be made clear. In the 
latter, interpolation is noccmally performed by constructing 
a polynomial fit of an adequate older to the given data. Often 
a linear interpolation suffices . A study of the process in 
the frequency dcmain shows that these methods are not the 
most suitable schemes for interpolation (1). Por the class 
of bandlimited signals, inteirpolation is shown to be 
essentially;, a low pass filtering process. Thus the objective 
in the desigi of an interpolator is to obtain the best 
approximations to. an ideal low pass filter. 
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A PPLICATIOITS OP DI&ITAL IMPIEPOMTIOII 

It is instructive at this stage to look at some of 
the applications of interpolation. Of the several applica- 
tions , we choose to highlight a few. 

(a) Yo coding ; 

In this method of data compression, information 
about slowly varying speech parameters obtained by the 
analysis of natural utterances is transmitted at low sampling 
rates. The regeneration of speech requires a higher sampling 
rate for synthesis ( 2 ), (3). This is achieved by interpola- 
tion. 

(b) Code format Conversion ; 

Por Delta Modulation (M) and Pulse Code Modulation 
(PCM), the required sampling rates are inherently different.^ 

A filter that effects a conversion from DM to PCM is 
basically a rate changing filter (4), ( 5 ), 

(c) Digital Realisation of Prequency Division Multiplexing 
(PlM)j 

In an efficient realisation of an PIM system employing 
digital filtering, complicated operations are performed at 
low sampling rates before translation to higher rates for 
grouping of daannels in the final stage ( 6 ), Here we have a 
need for interpolation at the transmitter and decimation 
at the receiver. 
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(d) T wo-dimensional Processing ; 

Sometimes the array of digital data in the sampled 
version of a picture may not be adequate for a good display. 
Interpolation helps put the picture into better evidence. 

On other occasions, for instance, in the data collected by 
satellite observation of the earth, the grid points in 
the input data represent a distorted perspective projection 
of the object of observation. Geometrical correction on the 
entire data requires excessive computations. Instead the 
correlation is performed on a sparse input grid and inter- 
polation is carried out on the corrected output grid, (32) 

(©) Efficient Pilter Realisations ; 

Interpolation techniques are now being put to wide 
use in yimproving ccmputational efficiencies of Digital 1 . 1’. 
Enters, The feature that is exploited here is that ■ . 

computations performed at a rate equal to twice the highest 
output frequency are adequate. If the sampling rate of the 
input is hi^, the data' cshh decimated and filtered before 
being restored to the original sampling rate by interpola- 
tion, A great deal of investigation has been carried out 
on multirate filtering (7)-(l2), 

REVIEW OE RESEARCH IN DIGITAL INTERPQIATiaU 

If we trace the developments in this field, we 
recognise the paper by Schafer and Rabiner in 1973 ( 1 ) 
as a landmark publication. It put the problem in the 
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righ-t sperspecti'^^ from the point of view of digital filtering, 
The authors suggested methods for obtaining change in sampling 
rates with rational ratios, A strong case was built up for 
the use of finite Impulse Response (flR) filters. It was 
soon recognised that when the initial and final rates differ 
widely, using a sequence of filters proves more efficient 
than using a single stage filter. Using this approach, 
a multistage procedure for rate reduction in the ratios of 
2:1 was proposed (5), (11), (13). Here the use of halfband 
filters considerably reduces the computations. Recently a 
design procedure for properly choosing a sequence of half 
band filters that meets specified fidelity criteria, has 
been presented (14). 

Without imposing conditions on the type of flR 
filters, a cascading procedure wherein the sampling 
rate reduction at each stage can be greater than 2 has 
been suggested in (15)^ The attempt is to optimally allo- 
cate the decimation to minimise the overall amount of 
computation. This approach has been extended to the general 
case of multistage, processing (16), The problem of design- 
ing a multistage interpolator for minimisation of storage 
rather than computation has also been examined (17), 

For interpolation, the use of Infinite Impulse 
Response (HR) filters in the conventional foim has been 
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shown to he less efficient than the use of PIR filters (l). 
But use of HR filters in phase-shifting networks for 
obtaining rate alterations in an efficient way has been 
proposed in (18), Another method utilizes a state variable 
description to produce realisations of sample rate reduction 
filters (10), 

Digital Interpolation has been analysed, in the time 
domain also and the analysis enables us to optimally design 
a filter for minimisation of error (l9),(20). In the 
optimisation, stationarity together with a knowledge of auto 
correlat|.on function is assumed. 

Another method of interpolation employs the Discrete 
Pourier Transform approadi (21), Por this method, interpola- 
tion kernels and expressions for error bounds have been 
derived, 

SCOPE OP THE PRESBUT WORK 

The present work is concerned mostly with the 
comparative evaluation of various interpolation schemes ^ 
using a representative set of signals and a performance 
measure based on mean square error, A new method of inter- 
polation by choosing filter coefficients on the basis of 
the received data has been tried. Without assuming a-priori 
knowledge of the autocorrelation function, an attempt is 
majde by using the Maximum Entropy Method (MEM) to estimate 
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the Autocorrelation function, This estimate in turn is 
then used to design an optimum filter. Also a direct 
operation on the time sequence data using an ICEM approach 
has been tried. 

After a preliminary discussion cn the interpre- 
tation of interpolation and interpolation schemes in 
Chapter 2, the suggested approach is elucidated in 
Chapter 3. A set of representative signals were used to 
test the efficacy: pf the interpolation techniques discussed 
in Chapters 2 and 3 and these results are presented in 
Chapter 4* Chapter 5 analyses the results an land makes 
recommendations regarding the application of the various 
schemes for interpolation. Listings of some of the 
programms lased in these studies are included in the 
appendix. 
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CHAPTER 2 

AH ANALYTICAL INqiQHT lETO THE 
IHTERPOIATIOH PROCESS 


It should he pointed out right at the beginning 
that we shall be concerned here only with linear inter- 
polation. In this chapter, we first take a brief look 
at the classical techniques of interpolation. \(e next 
develop a frequency domain interpiretation of digital 
interpolation and discuss design criteria for interpolators. 
A result that emerges from the discussion istthat EIR 
filters are particularly suited for interpolation. ¥e 
proceed to discuss methods for efficient implementation 
of ithe interpolation schemes. A novel IIR implementation 
described in (18), based on phase shifting networks, is 
also examined. 

In the latter half of the chapter, the problem is 
formulated in the time domain. The optimum filter for a 
class of signals with a known autocorrelation function 
is derived. 

2.1-1 Classical Internolation Techniques 

A vast amouint of literature is available on inter- 
polation techniques used in numerical analysis. ¥e 
restrict our attention to aspects that are relevant to 
digital interpolation. 
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Suppose that the values of a continuous function 
x'(t) are available at equispaced points in time. That is» 
a sequence x(n) is available to us where 

x(n) = x‘ (nT) (2.1.1) 

The requirement is to obtain a set of values at a different 
time spacing T* = T/I<» I integer. In other words we are 
required to estimate a sequence y(n) from the available data 
such that y(n) "'spproximataslxJ’tuTO. In numerical analysis » 
this obj ective is realised by constructing a polynomial of 
degree P which passes through C^+1) observations around the 
point of interpolation. The order ^ is chosen according to 
the accuracy desired. Clearly this interpolating polynomial 
is unique,. A convenient method of computation is to construct 
difference tables of various types like forward, backward and 
central differences (26). Based on these difference schemes, 
a host of int erpolation techniques like Lagrange, (Jauss , 
Stirling, Bessel and Everett have evolved. The form that 
is most amenable to digital interpolation is the Lagrange 
form. Por the pth interpolation, y(mL+#), between x(m) and 
x(m+l), let us use (Z+Z'+l) samples x(m-^K), x(m), ...» 

x(m+z'). Then 

y(mL+p) = £ L,‘.(mL+p) x(m^4J;) (2.1.2) 

k=-Z ^ 

where p = Q,. i, Ljj.(. ), the Lagrange interpola- 

ting polynomial' 

'■ ^ rf' ' ' "i? ' \ 'I . . 
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(.Xy.-X^) ... 

The polynomial as adapted to discrete time .• 

» IS given by 

(mL+p ) 

= (g-l+p/L) . , . (-(k-p+p/L) (.(k+l )+p/L) . . . 

(k+E)(k+K-l77.. (1) (^1) ... 


Clearly for p = 0 (2.1.3) 

Llc(ml) = 

which means that y(mL) = x(m). The input samples 

are thus preserved. The right hand side of ^ 

'>^•1.3) is in- 
dependent of m. Hence it is evident that ( o i \ 

2) can be 

implemented by a convolution operation using ^ pjjj fit 
whose coefficients are determined by (2.1,3) 

Domain Interpretation 

The analysis presented here is based on ^ 

^ ll ) . As in 

(2.1^1) we consider a sampled sequence x(n) 


x(n) = x’(nT) 

The fourier transform X' (-w) of x’(t ) is 


(2.2.1) 


(w) = J x'(t) e“3wt 


The z-Transform of the sequence x(n) is 

i' \ ‘ 




X(z)"= ^ x(n)2'"^ 

n=-(3D. 


(2.2.2) 


X(e3’«'), the Pourier transform of the sequence x(n) is 
related to X' (w), (22)» as follows: 

Z(e3''T) = l X X'U+ifS) (2.2.3) 

k=-oc> 

Consider x’(t) to be bandlimited, i.e. 

X’ (w) = 0 for iwj > B 

as shown in Figure 1. If T^7t/:p, then from (2-2.3) it is 
seen that 

XCe^'^) = I X' (w), -it/T < w (2-2.4) 

It is well known that the Sine functions defined as 

Sin(7ct/r) 

Sinc(t/T) = (2-2.5) 

(itt/T) 

can be used for interpolation. But such interpolation is 
not practicable as the Sine functions have an infinite 
duration with significant contribution even at points 
widely separated. Interpolation functions which decay more 
rapidly than the Sine functions, like 

fin(7it/T) 

I(t/f ) = — (2-2.6) 

Sinh(iit/!r) 

have been suggested (27)- For the purpose of implementation 
these functions have to be truncated. Instea.d of using such 

truncated f'unctions, it is more reasonable to look • for finit 
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Figure 2.1i (a) Fourier Transform of a continuotis time signal. 


("b) Fourier Transform of the sequence obtained 
by sampling with a period T ^ % /B * {v /T ) • 



Figure 2.2; Sampling rate increase (T’ = T/3, -w^ s= ic^, 

■' w* = ll/p* ). 

(a) Fourier Transfoim of sequences u^Cn) and x(n). 

(b) Fourier Transform of interpolated sequence y(n). 



duration functions. The search is facilitated hy the 
analysis that follows: 


2.1(a) Sample Rate Increase; 

Suppose it is desired to increase the sampling 
rate by an integer factor L. The new sampling period is 

T' = T/1 

Here we seek to fill in (L-1) samples between each two 
adjacent samples of the input data. A new sequence u^Cn) 
is defined as 

u^(n) = x(n/I), n = mL and m = 0, +1, +2» ... 

=0 , otherwise (2.2.7) 

The z-Transform of the sequence ^^(n) is 
oo 

O (z) = S;u.(n)z-“ 

n=-<30 

= X x(n)z"^^ 

= X(z^) (2.2.8) 

The Rourier transform of the sequence ■^.^(h) is 
Uo(e^'^') = Z(e^'^*^) 

= Ke^'^) (2.2.9) 

Uo(e^,'^ )'thus has a period 2ii/T rather than 2'n;/T^ as is 
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normally the case for signals sampled at a period 

T' . 

The desired sequence y(n) is such that 
y(n) = X* (nT ’ ) 

This is achieved .if we ansure that 

X' (w), -rt/T* ^ w ^ ti/T^ (2.2.10) 

' T’ • - 

The required filter sho\ild produce the sequence y(n) hy • 
operating on UQ(n). In the Transform domain 

Y(e^''^’) = H(e^'^’) U^Ce^^') 

= H(e^^*) xCe^"^) (2.2.11) 

That is, Y(e^’^^’ ) = i H(e^'^'^' ) X’(w), -it/T < w ^ n/T 

( 2 . 2 . 12 ) 

The implication of the above equations is illustrated in 
Figure 2.2. The requirement to remove the (1-1) images 
of (1 /T)X'(w) in U^(e^^*) which lie between w = n/T and * 
W = (2l-l)'rt/T and which are centred at w = 211/1, 4it/T, 

..., 2(l-l)it/T. A digital lowpass filter that rejects 
all frequencies in the range ir/T w ~ ic/T' meets the 
above requirements. From (2»2.10) and (2.2.12) we see 
that the passband gain must be 1. Clearly, the ideal inter- 
polator is 
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2Tt/T» 47c/T* 


Pigiare 2.3! Sampling rate reduction (T*=2T) 

(a) Fourier Transform of x(n), (b)Fourier 
Transform of y(n)vith. aliasing, 

(c) Fourier transfcm of y(n) obtained with 


a pre-cutoff . 



Figure 2.4: Sampling rate change by a factor 3/2. 

(a) Fourier transform of x(n) 

(b) Fourier Transform of interpolated inteimediate samples, 

(c) Foxirier Transform of final sequence. 
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=1 , iw/ < n/T (2.2.13) 

= 0 , ( 7 t/T)< M <(ix/T' ) 

2.1(b) Sample Rate Deer ease- Integer Factors: 

If the sampling rate is decreased by an integer 
factor M, then the new sampling period is T*=MT. The 
required sequence is 

y(n) = x(Mn) 

¥e therefore retain one out of every M samples with a "" 

spacing of pefit)d“T!. If the decimated samples are to 
uniquely determine x'(t), the Sampling Theorem requires 
that T ' n/B. If this requirement -is violated then 
aliasing errors occur. To avoid aliasing, lowpass 
filtering is performed before decimation. 

A sampling rate reduction with M = 2 is illustrated 
in Figure 2.3 for the two cases of T’ < n/B and T’ > n/B. 
Also shown i& the case where, in order to avoid aliasing, 
the sequence has first been passed through a lowpass filter* 
with a cuf-off frequency it/f * . 

2.1(c) Changing Sampling Rates by Ron-integer Ratios: 

We now consider changes in sampling rate by non- 
integer but rational factors (L/M). To achieve this, the 
above two techniques are combined by carrying out 
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interpolation by a factor L followed by decimation by a 
factor M. For the case M>L, if the new period T* is such 
that T' = MT/L '>71 /B? then aliasing has to be taken care 
of. Figure 2 *4 illustrates a sampling rate increase by a 
factor 3/2 • ,, 

2*3 Frequency Domain Design of DigLtal^ Int ernolation 
Filt ers 

It is clear from the frequency domain analysis 
that sampling rate alteration requires a lowpass filter. . 
An obvious .question that arises " ■ ' at this point 

is the type of the filter to be used. A basic considera- 
tion is the choice between FIR and HR filters. Other 
considerations rare regarding the specification of 
filters and the efficiency of implementation. These 
aspects are studied in this section. 

2 . 3 (a) The Case for FIR filters; 

The ideal interpolator should have a zero phase 
Or at most a linear phase. HR filters do not possess this 
property. On the other hand several techniques are known 
for realising linear phase FIR filters (23)» (24)* Thus 
FIR filters eliminate errors due to phase non-linearity, 
and also provide arbitrarily small values of error due 
to amplitude distortion. 

For conventional filters, HR implementation is 
inherently economical in terms of computational complexity. 



FIR filters generally requir© iftore computations for 
meeting the same specifications. These observations 
are true if the input and. output sampling rates are the 
same. But in interpolation, the nature of the problem differs. 
The output is obtained by filtering the sequence 
defined by (2.2.7) with a lowpass filter. Consider a zero- 
phase FIR filter with (2R+1) samples for increasing the 
sampling rate by a factor I. The choice of an odd ntimber 
of samples will be explained later. 

The output samples for the filter h(n) are 
n+Jff 

y(n)= SI U (J;)h(n-k) (2.3.1) 

k=n-N ° 


3ince only one out of every L samples of u^(ii) are 
possibly non-zero, the actual number of computations 
are reduced by a factor L* This feature is not exploited 
in conventional HR filters where the computations must 
be carried out at the output sampling rate. 

Again, in a decimation by a factor M, HR filters 
reqiiire computation at the input sampling rate as each 
incoming input sample updates the state variables of the 
filter. At the output (M-l) out of every M samples are 
discarded. On the’ confc rafy, considerable saving is 
achieved in an FIR filter by computing only one out of 
every M samples’. 
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For changing rates hy non-integer ratios, say, 

Jj/Ef we first increase the sampling rate by a factor L and 
later reduce it by a factor M. The features discussed 
above can be incorporated in each stage of the process. 

2.3(h) Impulse Response Constraints: 

One restriction on the Impiilse Response arises 
from the fact that the original samples should, be pre- 
served in the output. That is 

y(mL) = x(m) 

clearly the reqixirement is 

i 

h(o) = 1 (2.3.2) 

h(n) = 0 for n = +21,.... 

and }n| <ir 

Another restriction is tlmt the length of the 
Impulse response should be odd. It can be shown that for 
even number of samples, a linear phase FIR filter must have 
a delay of at least one-half sample (25)- Such a filter 
cannot preserve the samples of the original sequence, 
whereas a filter with an odd number of samples can* 

Yet another restriction comes in if we impose the 
condition that fo;^ the computation of each output sample, ±.e 
the number of input samples involved should be the same. 

For an, interpolation by a factor 1, the length of the 
filter response is chosen as 2EL-1. 
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2.5(c) Design Specifications: 

For interpolation by a factor 1, the approximation 
of the ideal interpolator consists of 

(i) specifying a passband within the range 0 s w < tc/T 

such that the gain over the band is close to imity. A margin 
of ^ allowed. 

(ii) specifying (L-l) stop-bands centered at 

w= (4it/T), ...» 2(lr-l)7t/P where the gain is 

close to zero within the margin 6^< 1. 


For cases where the bandwidth B is close to n/T, 
the specifications reduce to that for a lowpass filterf 
There are cases however of narrow band signals where a 
band stop filter of the same order achieves better 
results. This is because the filter coefficients can 
be suitably chosen to cause heavy attenuation only over 
small bands centered at 2m /I* 4m /I,..,. 


2»5(d) Efficient Implementation; 

For an FIR design, we have not yet discussed 
schemes which are computationally optimal. Row we see 

how some features in the interpolation process have been 

, ■ . ■ ■ 

exploited to reduce the number cf computations. 

Sv.,'!, 

Firstly thfe impulse response for a zero phase 
FIR filter is syimmetfic. That is 
h(n) = h(-n) 


(2.3.5) 
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This has inherent advantages for an interpolation 
hy a factor 2 as samples equidistant from the point of 
interpolation can he added first and then milLtiplied by 
the filter coefficient. Thus the number of multiplications 
is halved. This method suggests that if interpolation is 
carried out in a cascade of filters with 2:1 increase 
then considerable saving in computations is achieved. Such 
a procedure has been discussed in (5), (11) > and (13). A 
procedure which incorporates such stages has been presented 
in (l4)f where a set of nine half band filters are provided 
for combination to get the desired sampling rate change. 

A miiltistage procedure to minimise the computations 
using interpolation ratios which are not necessarily 2sl has 
bs#n described in (16)^ (Jiven the initialsampling rate, the 
overall interpolation ratio and the filter specifications, 
the overall Stm of the computation* rates at each stage 
is minimised. The flexible parameters for minimisation are 
the interpolation ratios I., i = 1,2,...,E as well as K, 
the nxmber of stages. Design curves and formulas for the 
implementation have also be eh presented in this paper. A 
similar approach has been used in (17) for minimising the 
storage rather than the n-umber of computations. 

2.4 Digital Interpolation bv a igplynhase NetWork 

In (18), a method of manipulating a lowpass filter 
response to make it amenable to efficient HR interpolation 
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has heen described. The authors hfifcire introduced a structure 
based bn phase-shifters and have referred tp it as a poly- 
phase network. Filtering is carried out by cancellation of 
the undesired signal components at the output of a smmation 
device connected to a set of phase shifters. The phase 
shifters are all pass filters and rotate the signal 
components in such a way that the components to be retained 
have the same phase at the sximmation point and add while 
those to be removed cancel on summation. If an intei^olation 
by a factor L is reqviired, that is^ the new sampling period 
is T' = T/I/» then the interpolator requires L phase 
shifters. 


Let H(z) be the transfer function a LPF which is 
obtained by conventional methods to meet the specifications* 
Suppose 


H(z) = 


A IT (z-z, ) 


K 

IT (z 


Pk> 


(1.4.1) 


The phase shifters are shown to be a cascade of a delay 
element and a filter which is a function of z^. The 
Polyphase Network transfer function is thus 

Hdzi= X 


(2*4.2) 
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To determine the coefficients of we equate 


(2*4«1) and (2.4*2). ¥e also note that 

l-l . I“2 I-l 

2 z + pj^ z + . • . + Pjj. 


z-p. 


L L 
z - Pk 


(2.4.5) 


Using this in (2*4.1) and equating with (2*4.2) 


KL 

a . z 

i=o 

k=l 


-n T 

s: 2 “ H (a^) 

n?=o 


(2.4.4) 


From a simple decomposition of the left hand side, it is, 
seen that 


=[ f (2"^)^-]/[n (1- pj a'^) (2.4.5) 

^ k=o r fc=l ^ 


where E* = Z for n = 0 

= E-1 for n > 0 

The fact that is a function of implies 
that the signal samples coming out of corresponding devices 
are weighted sians of input samples and previous output 
samples with N sample spacing. Moreover, from (2. 4.5. ) we 
see that the denominators are all identical. This is 
exploited in the implementation. Figure 2.5 illustrates 
the implementation for L = 2 and K = 3* 

The multiplication rate is (II)+E+l)f_/l where 
fg is the output sampling rate. This shows a considerahle 
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reduction on the 2K computations for the conventional 
HR implementation. Comparison with FIR filters should 
take into consideration the fact that filter coefficients 
for HR implementation for given frequency response 
specifications are fewer. The phase response howevery cannot 
be zero-phase for HR filters. 

2.5 Time Domain Analysis 

The problem of interpolation will now be formulated 
in the time domain using the approach of Oetken et al (19). 
The analysis has been carried out for noisy samples. 

Suppose we seek to increase the sampling rate of 
a signal by a factor L. The available r-sequence is x(n)f 
where 

x(n) = x' (nT ) 

The desired sequence ■u.(n) is 

u(n) = x'(nT')» T' = T/l (2.5.1) 

We recall the definition of uj^n) in (2i2#^7)' as 
u^(n) = x(n/L) = u(n), for n = ml' 

=0, , otherwise 

The problem is to design the interpolation filter 
h(n), which, operatic on the sequence UQ(n) gives an 
output y(n) to approximate u(n) according to some error 


criteria 
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We asstme that the filter impulse reponse is 
approximated hy = 2LK-1 samples where Z is an integer f 
K ^ 1. With this we obtain 
IK-1 

y(n) = SI h(i) u^(n-i) (2.5.2) 

i=lK+l ° 

Since u^Cn) = 0 for n^ml we get 

(i) For n = ml in (2*5.2) 

Z-1 

y(mL) = li(3l') u ((m-j)!) (2*5.3) 

j=-(K-l) ° 

(ii) For n = mX+ip; P = 1» 2> . • . ,1-1 

Z-1 

y(ml+p) = h(jL+p) u ((m-j)l) (2.5*4) 

j=-K 

Since we assumed hC-^Zl) = 0, (2.5*4) holds for p = 0. We 
note that for different values of p, different sample points 
of h(n) appear in (2.5.4). If we now consider the input 
sequences to be corrupted by noise, the contribution of noise 
to the output y(n) can be included in (2*5. 4.r,-) 

Z-1 

y(ml+p) = ST h(oL+p)[u ((m -3 )I) 4-w((m-o)L)] (2*5.5) 

3=-Z 

¥e now define I error sequences 

yp(ii) = y(n) - mCn) for n = ml+p 

= 0 , otherwise. 

f 

The sequence for total error is 
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I-l 

y(n) = i: y^(n) (2.5.6) 

p=0 ^ 

A reasonable performance index would be 

= ilyj! IlyJi ^ (2.5.7) 

p=o 

where ||.}| represents a norm. The norm can suitably be 
defined for various classes of signals. 

Since each error sequence y (n) depends on dis- 

!P 

joint subsets of h(n), the minimisation of i|fyjj can 

^ 2 

be done by minimising each of |fy j^i separately. 

Consider x(n) to be a sampled version of a stationary 
random process with zero-mean and autocorrelation function 
R(n). For this class of signals the performance measure 
is chosen to be 


® mi 

= E ([y(mL+p)-u(mI+p] 


(2.5.8) 


* b|[ 2 : h(jL+p) (u( (m-j )L) + w( (m- 3 )l)) - u(mlH‘p )3^ J 

’ ( 2 * 5 . 9 ) 

Assuming signal and noise to be uncorrelated 


3e* 

= 2[ r b(jI,+p)[lL ((3-1)1) + H ((3-1)1)] 

3h(il,+p)«A ^ ^ 

- B.^(il 4 -p)] (2.5.10) 
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where i = -E, . . . , 0 , . . . , Z-1 

R^(n) = B [u(m)u(m-n)] (2.5.11a) 

and R,(n) = B[w(m) w(m-n)] (2.5.1113) 

wf 

¥e define R(n) = R (n) + R (n). 

u w 

For minimm e^, from (2.5.10) we get 
E-l 

ST h(3l+p)R( (3-1)1) - R^(iL+p) 

3=-E 


For convenience we define 


R(o) 



Ed) 


Rd) 

R(o) 


. .. R((2K-1)L) 


R((2K-1)L) ... ..R(o), 





J 


h(-El+p) 


R^(-Kl+p) 

H = 

h(-KL+l+p) 

• 

« 

• 


R^(-KL+l+p ) 

• 

« 


h((K-l)l+p) 

i 

[ 

R((E-l)l+p). 

„ J 


(2.5.12) 

f hen we have 

S " i * (2.5.13) 

The computation of filter coefficients is thus 
seen to req,uire the inversion of a 2Ex2E matrix and 1 
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multiplications iDy the vector The symmetry of the 

matrix can he exploited hy adding the first and last rows, 
the second and second last row and so on. This operation 
allows us to break the system of equations into two sets 
which require the inversion of two KxK matrices. 

For a noiseless signal the choice of h(mL) is the 
same as in (2.3«2). 

The optimum filter thus designed serves as a standard 
for Judging the performance of the other filters. Expres- 
sions for error are derived in Appendix I. In the 
evaluation of the various schemes, the errors will he 
studied in the context of the optimum interpolator. 



CHAPTEE 3 


DESIGN OP DIGITAL IM) EBP PLAT 0R3 USING- 
AN ESTIMATE OE THE SPECTRIM 

In Chapter 2» we had examined approaches to digital 
interpolation of a signal about which no additional infor-- 
mation is available. We also saw that with a knowledge of 
the autocorrelation function, an optimum filter for minimum 
mean square error can be designed. In this chapter, we 
investigate the possibility of improving on the digital 
interpolation of a signal about which no apriori information 
is available. The approach suggested is to estimate the 
spectrum and design an optimum filter using this information. 
Essentially the idea is to compute the autocorrelation 
f^unction from an estimate of the spectrum and design an 
interpolator for minimum mean sqxiare error. 

The above proced'ore will be discussed in detail 
in Section 3* 1* Section 3»2 we shall look at the 

optimisation procedxire in the frequency domain. This helps 
us to get a better insight into the optimisation. The 
final results are the samer as those obtained in Section 2.5* 
In Section 3.3, we briefly mention the methods of spectrum 
estimation with the focus on the Maximum Entropy Method 
(MEM). Finally we study MEM interpolation by direct operation 
on the data as suggested in (28). 
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5«1 Procedure for Obtaining the Approximate Optimiam 
Interpolator 

Given the sampled version x(n) of a continuous 
signal » with a sample spacing T, we again consider the 
prohlem of digital interpolation hy a factor 1. In order 
to design an optimum interpolator as in 2. 5> values of 
autocorrelation function at a lag T' = T/1 and its 
m-ultiples are required. Clearly by a direct operation 
on the data, we can, at best, estimate the autocorrelation 
function at a lag of T and its multiples. To overcome this# 
we first estimate the spectrum and evaluate the autocorrela- 
tion ftmction at the required resolution. As usual, we assume 
that the signal is bandlimited. 

Using the input data we first calculate the spectrum. 

¥e use the MEM estimate in view of its recommendation for 
being non-comittal about data outside the record length 
and its adaptability to slowly varying statistics. Having 
obtained an estimate of the spectrm, we can find the auto- 
correlation f-unction at any desired lag. We can use the PPT 
algorithm for computing ACP. The next step is to employ the 
procedure described in Section 2*5 to obtain the optimtm 
interpolator. 

How the spectrum estimate helps in bringing about 
an improvement in performance will be clear from the next 


section 
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3 • 2 Frequency Domain Interpretation of OptimlBation 

In Chapter 2 » we had seen that interpolation of a 
signal hy a factor L implies that in the frequency domain, 
we must remove the images of X'(w), the Fourier transform of 
the continuous signal x'(t), centered at w = (2'n:/T)» 

(4'n;/T), ... 2(L-l)7i/T. Refer to Figure 3*1» I'or the 
purpose of interpolation, these images can he viewed as 
noise. The problem is thus analogous to the classical 
problem of designing a filter for minimiAm mean square error. 
In this case the noise and signal occupy different regions 
of the frequency band. 

For optimisation we require the definition of 
an error measure. This could be achieved by writing the 
eqtii valent frequency domain expression for the error in 
2.5» We shall obtain the same expression in the zifollowing 
manner. 

x(n) is the available sequence and nQ(n) is the 
sequence obtained by inserting (L-1) zeroes between each pair 
of adjacent samples of x(n). The autocorrelation function 
sequence for x(n) and UQ(n) is 

R(n) = l(x(R) xCl+n)) (3*2.1) 

R^(n) = 1(u^(K) u^(Z+n)) 

Clearly Ro(n) is the sequence obtained by inserting (L-1 ) 
zeroes between each adjacent pair of samples of R(n). So 
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Eigtire 3*1 Interpolation by a factor L =s 3 , 
Hj : I(i.eal Interpolator 
H s Approximation to 
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the Fourier transforms of these two sequences are the 
same. ' 

= 3 ( 6 ^”^) ( 3 . 2 . 2 ) 

Let Hj(e^^ ) be the ideal interpolator and H(e^^ ), the 
approximated interpolator. The error power in the frequency 
range w and w4dw due to the difference between Hj and H 
is (but for a constant factor) 

dpg = SQ(e^'^'’)/Hj(e^'^*)-H(e^"^')/2 dw (3.2.3) 

Summing the error power we can express our error meas'ure 
EM as 

EM = 5^*^. S(e^^^)/Hj(e^^* )-H(e^^' )/^ dw (3*2.4) 
o 

¥e have seen that 

H^(gdwT,») ^ ^ w Sit/T 

B= 0 for It /T < {w| < ti/T ’ ( 3 * 2 . 5 ) 

Therefore, 

M = J ^ S(e^'^)/L - H(e^^')/^ dw 
o 

+ r^' /H(e5'^')/2 s(ej“^)dw (3.2.6) 

‘'ii/T 

¥e thus see fhom (3.2.6) that the spectrum serves as a 

i ig. > * ’ ■ ' ^ 

weighting factor for the error in approximating an ideal 
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interpolator. The optimum filter would be the one which 
has a response close to the ideal interpolator over 
frequency ranges where the signal power is significant. 
Now, for linear phase, 

h(n) = h(-“n) 


Using the notations as in Chapter 2» 

H(ej^') = t: n 2 CosC (iL+p)wT ' ] (3-2.7) 

i=o p=o - 

where a^ = ih(o) 

and a^ = h(iL+p) p= 0,1,... ,1-1 (3-2.8) 


,7t/LT 


K-1 L-1 


EM = I ' -[L-2 SZ IZ % Cos[ {iL+p)»T']]^ S(e^'^^)dw 


1=0 p=o 


%/T' K-1 L-1 


' f ' C2--iI ^ a* GOB((iL+p)wT' )]^ S(e^^)dw 
^ it/LT’ ^=0 p=o ^ ' 

(3-2.9) 


d 


9 ^^ 


# = 0 for minimum error. Prom this condition we get 


®1 

E-1 Ir-1 . n/T' ^ 

2 51 21 / Cos[ (i.L+p^)wT’3S(e^^^)dw 

i=o p=o ^ o • - 


. 1.0 


lij Cos( (i^^L+p^wT )S(e^^)dw 


K -1 1-1 


n/H 


^ 21 ' Cos((i-i^L + p^j_)wT' )S (e^^^)dw 


( 3 - 2 . 10 ) 


1=0 p=o 


o 

%/T^ 


+ jf ' .COs((i+i^L + p+P 2 ^)wT* ) S(e^^^)dw] 
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%/m 


L f Cos( (i3_L+p^)wT’ )S(e^^^)dw (3.2.11) 


is periodic with, period 271/1. Hence the first 
integral on the left hand side is non-zero for (p-p^’) = 0 , 
And the second integral is non-zero for (p+p^) = cil. We 

therefore write (5. 2.11) as 

K-1 . .7t/T‘ . T, 

^ %1 f •Gos((i-i3_)wIT' )S(e^^^) dw 

i=0 - n 


/ K-1 . 

+ H ai / Gos( (i+i3^+l)wLT’ )S(e^^)dw 

i=o o 


%/!! 


j Cos( (i^l+Pl)wT ' )S(e^'^^)dw 


(5.2.12) 


where p^ 0< 


(3. 2. 12) is identical to the results obtained in Section 2.5. 
For p=0, this also yields 


and 



for i = 1,2, 


,K-1 


(3.2.13) 


let us apply the optimisation to a single tone input with 
a random phase 6 


x’(t) = Gos(wQt+6) 

E(t) = Gos(w^i:) 

S(e^^^) = ^ 6(w-Wq) 0< w <7t/T (3.2.14) 

Let 1=2, E=2. For p=o, the coefficients a^ are given by 
(3.2.13). Foi* p=l» we have from (3.2.10) 



( 3 . 2 . 15 ) 


TL 2a^ Gos((2i+l)wQT» ) = 1 

i=o 

We see that one of the two coefficients, a^ or a^ can he 
a.rhitrarily fixed. This is apparent as we are require 
to choose the coefficients so that the gain at the point 
w=Wq is 2 and is immaterial elsewhere. Choosing a^ = 0, 
we get 

a° = l/(20os (w^T ' ) ) 

H(e^'^^V= 1+ — i — CoswT' (3.2.16) 

Oos(WqT*) 

Clearly at w=w^ the gain is 2 as reqmred. For narrow- 
band signals it is thus sufficient that the gain he L over 
a narrow hand of frequencies. ¥e ca,n also conceive of an 
adaptive interpolator for signals with a slowly varying 
'spectrum*. 

Another point that is clear from the frequency 
domain interpretation is that for 1=2, the optimisation 
simplifies considerably. For this case the pass-hand 
and stop hand are symmetric and equation (3.2.12) assumes 
a simple form with = L-p^. 

¥e now briefly divsouss the methods of estimating 


the spectrum 



37 


3*3 Spectr’um Estimation 


Among the conventional methods of spectrm estima- 
tion, Blackman-Ttikey and Periodogram methods are popular. 
The Blackman- Tukey approach first estimates the auto- 
correlation function and then expresses the spectrum as 
a Pourier transform of this sequence. Auto covariance 
estimates outside this range are assimied zero. In 
Periodogram technique, the following estimate is used 


S» (w) 


1 


N 

X x(n) 

rt=l 


2 


(3.3.1) 


Here the assumption of a periodic extinsion of data is 
implicit. These methods suffer from several disadvantages 
apart from the assumptions in extending the data. 


in recent years, the Maximum Entropy Method (MEM) 
of spectrum estimation has been recommended for its 
several advantages. We shall not go into the theoretical 
details of this method. Of the several procedures, we 
shall use the popular Burg’s recursive' algorithm (29) for 
fitting in Autoregressive Model to the input data. This 
method minimises the error power in passing the input 
series in the forward and backward direction through the 
filter. The optimum order of the Autoregressive Model 
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is determined by Akaike’s Pinal Prediction Error 
criteria (50)* Por a detailed study of methods of 
spectrvim estimation with emphasis on. MEM estimatiion, 
the reader is referred to the thesis (3l)» This thesis 
contains useful detailed discussions^ programs and 
extensive bibliography on the subject. The flow chert 
for the procedure is given in the Appendix II. The 
details will be clear from the routine SPEMIM in the 
Appendix. 

3*4 Direct MPM Int ernolation 

This method is suggested in (28). The basic 
idea is to determine the Pourier Transform of the 
input signal using the Maximum Entropy Method. Prom a 
knowledge of the Transform, the signal samples can be 
evaluated at any instant of time. 

By a classical statistical mechanics argument, the 
author has related entropy of the sequence to the Pourier 
Transform by the following equation 

H ln[ X(w)2^(w)]dw (3*4.2) 

-B 

This is maximised subject to the constraints 

f X(w) df = x(n) 

where w = 2nf 


(5.4.3) 
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and n= -M, 0, M. The solution to the 

above is shown to be 

X(w) = [ I: (3.4.4) 

n=-M ^ 



. ,M 

,M 



(3.4.6) 
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Once the coefficients are determined, the Fourier 
Transform can be computed at any desired resolution using 
(5*4.4)« The interpolation can be carried out using 
(3*4*3)* Clearly in view of the matrix inversion as required 
in (3*4*6), this method has severe limitations with regard 
to data length. Also, the number of computations per 
interpolated sample is extremely high compared to 
conventional interpolation. 

A routine MEMIET for performing interpolation 
using this technique is included in the Appendix. 
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SHATTER 4 

RESULTS OF STUDIES USIM A REpRESEOTATIVE 

This chapter presents the results of application 
of the interpolation techniques discussed in Chapters 2.. 
and 3 to a representative set of signals. In Section 4»1 
we give a description of the various signals employed for 
staying the performance of tho interpolation techniques 
and the filters used for inteipolation. In Section 4*2 
we list the results of studies on deteiministic signals 
and in Section 4*3 on random signals. 

4.1 Description of the Signals and Filters used for 
the Study 

For a meaningful comparative stiJdy, it is necessary 
that the conditions and restrictions oh the various digital 
interpolators he identical. In this study the initial 
sampling frequency is chosen as 1 Hz. The signal is consi- 
dered handlimited to 0.5 Hz. The factors hy which the 
sampling rate is sought to be increased are chosen to be 
L=2 and 1^=4. Four neighbouring points are used for inter- 
polation of each sample. Thus the FIR filters used have an 
impulse response of length 7 and 15 for the cases Ii-2 and 
L=4 respectiv^y. The input record length was chosen as 128, 
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4- 1(a) Signals used: * 

Only real input series will be considered. The 
performance will be evaluated on the basis of the mean 
square error. For obtaining the error, the samples of the 
continuous signal at the points of interpolation i^hotild be 
available. The approach in the study has been to first 
obtain a series at the desired final sampling rate. This 
series is then decimated to get the input data for inter- 
polation. The error sequence is thus the difference 
between the interpolated sequence and the undecimatai 
input sequence. 

Both deterministic and random signals have been 
used in the st\idy. Single tone inputs, sums of sinusoids 
and Sine functions are the deterministic signals used. The 
random signals were generated by processing a White noise 
series, generpted by the method of multiplicative noise 
(routine WGNOIZ). The white noise series was passed 
through a single first order lo^ass filter or a cascade of 
these filters. Abother series was obtained by simmaing up . 
24 consecutive signal samples in the white series. This 
series has a triangular autocorrelation function and is 
useful for its analytical tract ability. 

24 

x(k) = 21n(k+i-24) (4.1.1) 

where n(k) are white noise samples. 
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The autocorrelation function is 


Il(n) = I!(x(k)x(k+n) ) 


(4.1.2) 


R(n) = 0^(24 ■“ n ) for n < 24 


(4.1.3) 


= 0 


for n >24 


where E(n(K)n(Z+n))= (4. 1.4) 

4.1(b) niters used: 

The notations in the brackets will be used to 
refer to the filters. 

(i) Truncated Sine fxmetion (TSIRC): In this case the 
finite duration impulse response is a truncated sine 
function extending over 4 seconds, i.e, four sampling 
periods of the input signal. 

(ii) Cubic Spline Response (CSPLIN); A finite duration 
cubic spline curve used for approximating a sine function 
is mentioned in (32). The conditions used in the approxi- 
mation are a duration of four sample periods, slope and 
value continuity, symmetry about the time origin, and a unit 
area under the curve. The resulting cubic convolution 
polynomial is 

, 1-2 \t|^. + It} * , |t| < 1 

h(t) = 4-8jt|+ 5|tl* - |tp , 1 < |t| < 2 

■- “ r'*( '• >c >* 2 


(4.1.5) 



(iii) I-agrange Interpolator (IGRIJ),: a four point 
Lagraiige interpolator h.as beep used in tlie study* The 
impulse response samples for cases (i ) , (ii) and (iii) 

9-^0 listed in Table I. The coefficients for each inter- 
polated sample are given ro'w^wise* The impulse response 
and the corresponding frequency response are shown in 
Figure 

(iv) -Optimum Interpolators (OPTIM): We have used four 

point optimum interpolators derived from either a knowledge 
of the autocorrelation function (OPT I) or by the use of 
MM spectrm estimate as discussed in Chapter 3 . (OPTMM). 

(v) Polyphase Network (PPNET): IIP interpolators using 
phase-shifters were discussed in Section 2.4. This 
technique has been employed Ifbr d oubling the sampling 
rate with a third order Butterworth filter. The details 
of the filter coefficients are provided in Table II. 

See Figure 2*5 and Section 2.4 for details. 

(vi) Direct MM Interpolation (MMINT): The method outlined 
in Section 3«4» has been applied for interpolating some 
signals* In this case the record length was small* 

4 * ^ Results for Deterministic Signals 

The criteria used in evaluating the performance 
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TARLF J. I'^PULSE RESPONSE COEFFICIENTS 


(1) TRUNCATED SINC FUNCTION INTERPOLATOR 

G,OnOOuOOO 0«00000000 1.00000000 0.00000000 

“0.12T61666 0.30010548 0.90031632 -3.18006324 

-0.21220660 0,63661980 0.63661980 -0.21220660 

-0,18006324 0.90031632 0.30010548 -0.12861666 

(2) CUBIC SPLINE INTERPOLATOR 

o.ooauoooo o.ocoooooo i.oooooooo a.oooooooo 

-0.04687500 0,29687500 0.89062500 -0.14062500 

-0.12500000 0,62500000 0.62500000 -0.12500000 

-0.14062500 0.89062500 0.29687500 -3.46875000 

13) LAGRANGE INTERPOLATOR 

‘ : ■■ ■ ■: ■ ■ ' ■ ■* ■ -feSI*' ■ V ■■ ■■ *'■ - ‘ 

0.00000000 0.00000000 1.00000000 ^ 0.00000000 

-0,03906250 0,27343750 0.82031250 -0.05468750 

-0,06250000 0«56250000 ,0.56250000 -0.06250000 

-0.05468750 0.82031250 0.27343750 -0.03906250 




(a) Truncated Sine 
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Xu^(n) 

(S/N) = — (4.2.1) 

;^[y<n)-u(n)] 

where y(n) is the interpolated series and u(n) is the 
undecimated input sequence. The summation is over all 
interpolated samples. 


4. 2*1 Single tone input: 





‘i^ate by L=4 in 
4.^c>n was i|sed twice ■ 
age. The resialts are 


liliR 


listed in Table III are results for a sampled 


sinusoid of frequency O.J0^5 Hz. 'Using the approach of 

..i ^ . 'I, i .. M '>' / 

Section 5.2» the theoretically designed optimum filter for 

V \,a’' ■' ,/ « ^ jXi 

the known frequency was used (1=4). , 


TABLE 

III: Single tone ir 

iput c(f=?t). 0 Gi^ Hz) : 

,r- ' !•' ”, r 


Eilt er 

"■"T - — TT— 

'* |S X > -f ^ 

I'i 1 !' ’• 7 ', ■■ ' •' , ' 'j, '■ 

(S/H) 

t V . 

1 

2 

TSINC ' ' 

CSPLIH 

if w 1 

4 hv'--'fe ' 


3 





4 

.4. . . , ’ T 

OP TIM >■ 
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Table IV: Single tone input 
(f = Hz) 

Filler. S/N 

1. Single stage IGRNJ 199 *7 <3^1 

2. Two stage IGRHJ 201.8 db 


Interpolation Of a single tone* input by a factor 1=2 
was carried out by using a polyphase network. For an 
input tone frequency of 0.03125 Hz, a phase lag of 0.025955 
radians was observed. The error was computed after elimi- 


nating the phase lag. The result'^^ is given in Table V. 

1 » fr**' 1 » I >tf<. ' *Vi S* ‘ *' ffUl i«V*t i,# W 1 I > <*T *> •< * 



Table V: Single'- tone input 


,v ,:.^52*5 ^ 

. ■' pt'T 

Direct MEM Interpolation is not applicabl 
The solution matrix for a periodic signal 


Hz with equal amplitudes 


e 
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If the frequencies are known, an approach similar to that 
given in 3*2 can he used in designing the interpolator 
for exact interpolation. In the present study, the 
spectrum was estimated hy using the MEM approach with a 
maximum AR model order 10. The spectrum estimate is 
shown in Figixre 4.2. 

The theoretically computed impulse response 

« 

coefficients for exact interpolation for a factor L=2 
and those obtained by MEM spectrum estimate are given 
in Table VI. 


Table VI; Filter Coefficients for Intermediate 
Points with L =2. 


Filter 

1 -1 
ai - 

a^ = 

1. Theoretical 

Opt im\mi 

0. 6199144 

-0.1920603 

2. OPTMEM 

0.6208004 

-0.1921782 


A comparative st\3dy of the performance with a 
Lagrange interpolator can be made from the following 
table. 
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ligiire 4 *2: MEM ^ectriim for a sum of Sinusoids* 



li-gure 4*3J Frequency response of t lie optimum 

filter for a sum of sinusoids obtained 
by an MM spectrum estimate. 

= 0*125 Hz and f^ = 0*333 Hz). 
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Table VII: Sum of two sinusoids (Ip=2) 

(f^ = 0.125 Hz and fg = 0.355 Hz) 

Filter (S/N) 

1^ OPTMEM 109.6 db 

2^ ..IfGRlTJ 38.3 db 


The filter response is shown in Figure 4.3. The 
filter gain at f^^ and f^ is 2.0000013 and 2.0051568 
respectively against the required value of 2.0, 

Another study was made using a sum of five sinuSoidfil 
closely spaced in frequency to get a narrowband signal. 

2 

x(n) = 51 Sin(2Kn(0.145 + o.Ol i)) (4.2»2) 

i=-2 

The restilts are listed in Table VIII. 

TableVIII Results for smi of 5 sinusoids 
(Eqn. 4.2.2) 


Filt er 

1 TSINC 

2 , CSPIIir . , 

3 LGRRJ. 

4 C^TIfflf" 

a , . , yp.iHpi<pi . 


(g/H) 

67. 0 db 
61.4 db 
9Q,7 db 
145.91 db 

CENIHaL urn Am 

Acc. No. A 
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4.2(c) Sine Function: 

The sampled sine function serves as a convenient 
deterministic function for studying direct MEM inter- 
polation. A 9-point input was chosen to cover the 
central lobe and the first sidelobe on either side. The 
results are listed below. 

Table IX; Interpolation of a sampled Sine function (1=2). 



Filter 

(S/N) 

1 

TSINC 

48.4 db 

2 

CSPIIN 

64.1 db 

3 

IGJRNJ 

112 . 9 db 

4 

MEMINT 

132.0 db 


4.5 Handom Signal Interpolation 

Weyfijrtjtt. consider the signal described by eqn. ( 4 .I.I) 
for which the autocorrelation fvinction is triangular. The 
optimum interpolator for this case is the linear inter- 
polator. The results for this case are listed in TaHe X. 


Table X; Interpolation for Signal with a 

Triangular Autocorrelation function. 


. 1 . ii..., ...I... 1 p.— . 1 1 

Filter 

(S/H) 

TSIIO " ' 

15.66 db 

C^lIN 

17.84« dt- 


18.99 db 

IINEAR ’ IHTEBEOIAT OR 

19,18 db 

. ' XOPTMBM 

17.95 db 
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. Plgure 4.4s ^ectr'oii for A = 0.9» B = 0.1. 


, r t ^ ‘ 
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Next we consider a white noise series passed 
throiigh. a first order filter of the foim 

x(k) = A x(k--l) 4 - Bn(k) (4.3.I) 

Second and third order filters were realised as cascades 
of identical first order filters. The normalised spectrim, 
with A = 0.9 and B = 0.1, for first, second and third 
order filters is shown in Figure 4.4. 

¥e shall now list the results for white noise 
series i Both uniformly distributed and (Jaussian, passed 
throvigh second and third order filters. ¥e first consider 
white noise passed through the following filters . (Tabl e XlfXIl). 


Table XI: SeGond Order Filter 
(A*0.95. B = 0 . 05 ). 



Filter 

(S/N) 

1 . 

TSINC 

41«37 db 

2 . 

CSPLIN 

65,75 db 

3 . 

LGRNJ 

69.92 db 

4. 

optmm"^ 

70.27 db 

5 . 

±4 

OPTMM^ 

(69.68}db 


OPTMEM^ : This is the optimt®, interpolator derived from 
the MM estimate computed at ft low frequency resolution 
(1/4 of that used in 4)* 

The restil:i® w!3ich will, be listed were 
obtained for Gaussian white noise passed through second 
and third order filters. (Table XIII - XVI). 




Table XII: TMixi Order Filter 
(A = 0.9, B = 0.1) 



Filter 

(S/N) 

1. 

TSIHC 

(39.70 db 

2. 

GSPIIN 

67.36 db 

3. 

LGilNJ 

84.19 db 

4. 

OPTMEM 

85.34 db 


Table XIII: (A = 0.8, B = 0.15) 


Filter 

Second order 
(S/N) 

Third Order 
(S/N) 

1. TSINC 

29.63 db 

36.81 db 

2. CSPIIN 

31.31 db 

51.57 db 

3. LGRNJ 

32.33 db 

58.92 db 

4. OPTMM 

' 32. 20 db 

59.46 db 

Table XIV: 

(a = 0.85, B = 

0.1) 

Filter 

Second order 
(S/N) 

Third order 
(S/N) 

1. TsiNc'"'”'' 

"'"'*''"■54.55" db ' ■ 

42.50 db 

2. GSPLli^ 

?-■: ■5;,r; i. 

43.68 db 

66.00 db 

3. LGRITJ 

44.32 db 

lul ^ '■ ' ' 

74.84 db 

4, OPTJSEM - 

44. 65 db 

75.60 db 
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In the following table, we list results for the 
cavse where the sample realisations used for computing the 
spectrum and for interpolation are different. That is, 
the white noise realisations that were passed through the 
third order filter were different. 


Table XV: (A = 0.9, B = 0.1) 



Filter 

Sec ond ord er 
(S/M) 

Third order 
(S/M) 

1. 

TSINC 

41.51 db 

41.12 db 

2. 

OSPLIN 

62.87 db 

77.82 db 

3. 

LGRNJ 

65 . 90 db 

86.27 db 

4. 

OPTMM 

66.56 db 

86.51 db 


Two seven point series were used for direct MM interpolation. 
In one case the input samples were preserved in the output. 

In the second, the output was not consistent with the input. 
Por the first case, the result was 

Table XVI: Second Order Filter 


(A = 0.95, B = 0.05) 


, — — 

Filter 

(S/M) 

MMINT 

55.28 db 

■ ’ " '" iGRNJ ' “ 

59.24 db 









CHAPTER 5 
C0HCLUSI0H3 
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In this chapter we analyse the resiILts presented 
in the last chapter. The results for the deterministic 
signals will he discussed only to the extent that they 
shed light on the processing of random signals. On the 
basis of the performance, some comments and recommenda- 
tions regarding the various interpolators are made. 

In our studies, only a four point interpolation 
has been used. This condition imposes a restriction on 
the approximation to the lowpass filter that can be 
obtained. In order to make a comparative study, we choose 
the Lagrange Interpolator as a standard for evaluation. 

The Lagrange interpolator is an approximation to the ideal 
interpolator with a relatively flat response at the origin. 
In the absence of any knowledge about the signal, it would 
seem that the Lagrange interpolator would be the best 
solution for interpolatioi^. This observation is confirm®3 
by studies on signals with frequency content veiy low 
compared to half the sampling frequency. 

In studying the interpolation for a single tone 
input, the optimum interpolator was derived assuming a 
knowledge of the frequency. Exact interpolation is possible 
here as the optimum filteris-essentislly a solution to a 
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set of equations- with two unknowns, the phase and the 
amplitude. A comparison of the other filters with this 
optimum is not meaningful. But the results show that the 
Lagrange interpolator performs well at the frequency 
f =0.00625 Hz. (Table III). The TSIHC and CSPLIN perfor- 
mance is far inferior. This is to be expected in view of 
their frequency response close to the origin. 

In the interpolation of a sum of two sinusoids 
the frequencies chosen were high (1/4 and 2/3 of half 
the sampling frequency). Here instead of assuming any 
knowledge of the input, we used the suggested method of 
estimating the MEM spectrum and designing the interpolator. 
The filter coefficients computed are quite close to the 
theoretical optimum (Table 71). The results in Table 7II 
show the wide difference in the performance between OPTMM 
and. IGrRHJ interpolators. Thus for signals with a large 
high frequency content, the spectrum estimation enables us 
to design interpolators for far better performance. The 
implications for interpolating a periodic signal are also 
clear. A gain of L at the various harmonics is adequate for 
exact interpolation. 

Next we look at the results in Table 7III for a 
sm of five closely spaced sinusoids. Here too the 
OPTMM interpolator is far superior to LG-RNJ. The gains 
Of these two interpolators^, at the mid-frequency cf the 
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sinusoids were 4.000839 and 3*963334. Thus if the 
estimate indicates a narrow band of frequencies, the 
interpolator can be designed with considerable saving in 
filter length to meet the ideal interpolator requirements. 

If we desire to increase the sampling rate of a m-ultiplexed 
signal, say of five composite signals, then the optimum 
interpolator would be required to approximate the ideal 
interpolator in the region around the five carrier 
frequencies with attenuation at the unwanted images. 

From Table IV we see that the two stage Lagrange 
Interpol at or ■-performs slightly better than a single stage 
filter. The frequency response of the cascaded system at the 
input frequency is closer to the ideal interpolator than 
the single stage system. This shows that a cascaded 
system can be j'udiciously designed so that the final 
gain in the frequency band of interest is close to L, 
the ■;'-rit :.interpolati on ratio.:, • This achieves a consi- 
derable reduction in computations. 

Table V shows that polyphase inpleiientation 
performs well but for phase delays In cases where 
phase distort ion,, is not important, e.g. speech signals, 
polyphase implementation can be used. An advantage in 
computiationw, results if the interpolation ratio is high. 

In studying thfe results for randaa signals* we 
shall first check how close to the theoretically expected 
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nininun error perfomance are the results actually obtained. 
Using the method of error computation in Appendix I, the 
(S/U) ratio for the optimum., in this case, a linear inter- 
polator was found to be 19.20. The value obtained from 
the implementation is 19.18. SPor this signal the OPTMEM 
performance is inferior to the Lagrange interpolator. This 
is because the signal spectrum has considerable high 
frequency content and this shows up as aliasing in the 
spectrum estimate. The optimum interpolator design thus 
uses an incorrect weighting in computing the filter 
coefficients. 

Now we consider the signals where white noise is 
passed throxigh a second and third order filters. In all 
cases except one, the OPTMIM performance is superior to 
lagnenge interpolation. This has to be seen in the light 
of the fact that that spectrxim is significant only in 
regions close to the origin. In this region, the 
Lagrange Interpolation approximates the ideal interpolator 
very closely. The exception is in the case of the second 
order filter (Table ZIII), where there is some aliasing 
in view of the choice of filter coefficients. 

It is important that the spectrum be evaluated at 
a . sufficiently high resolution as a computation using a 
low resolution spectrum may introduce errors in filter 
coefficients. This may be seen frcmi Table XI. 
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If the spectra have a considerable high-frequency- 
content then the perfomance could be improved greatly by 
using the spectral estimate. If the signal statistics are 
slowly varying, the spectrum can be evaluated periodically 
and interpolation coefficients updated accordingly. 

So-faraasMEMIOT is considered , its application required 
that to get a close estimate of the spectrum all frequencies 
should show up in the record. This would demand a large length 
Of the record. On the other hand, the solution matrix 
becomes increasingly difficult to solve with larger lengths. 
Por symmetric data the performance is found to be consistent 
but in other cases some inconsistencies showed up. Prom 
the point Of view of computational simplicity, this method 
is hardly competitive. 

We find that the Truncated Sine and cubic spline 
interpolators invariably showed a poor perfomance. Another 
result that emerges is the consistently good performance 
of lagrange interpolator. But the OPTMM interpolator perfoi^- 
mance shows that an improved design is possible with an 
estimate of the spectrum. This method is particularly 
usef\a for stationary random signals with large record 

length. This method can be advantageously used for narrow 
band signals. 



APPBKDIZ I 


62 


EXPRES^IOF FOR IRTFRPOLATIOR ERROR FOR STATIONARY 

SERIES 

Tlie error measure is given by (see eqn.2.5.9) 

a r 

®p ~ ^ R w((m-j)L) - u(ml4-p)]^^ 

(A.l.l) 

Ass-uming that the noise is uncorrelated with the signal 
2 K-1 K-1 

®p = j^+P)ii(iI'+p)[E^( (i-j )1) + 


+ \((d-i)P)] 


For minimum error 


E+1 

+ \(o)- 2[ .5Y h(jI+p)R^(jL4-p)J 
D — E 

( A • 1 • 2 ) 


3h(iL+p) 


This gives ' 


Z-1 


^ ^(3^+p)E( (D’'i)P) = R^(iE+p) 

(A.1.3) 

where 


. E(n) = R^(n) + R^(n) 

(A.1.4) 


Therefore, 


®p min ^(3^+p)Ry^(ji-tp) 
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Let a.^ = 

" \((J-i)L)/Il^(0) 
^ip = \(^^+i>)/\io) 




L-1 K-1 

21 [1 + X 


K-l 

21 . 


p=0 i=-Z j=:-K 


li(iL+p)h(jI+p)[a. .+§ . .] 


Z-l 

-2 ll(3L+p)y 1^-1 

0=--K 



, L-1 Z-l 

i pS j^-K ’^minfJl+PJYipr^ 




APPMDII II 


COyipUTE , 

J = 0,1, . . . , 


CCMPDTE EPB(M) | 


’B <MIN 


AMI]I(J)=a(J) 
MlfcpPE(M) ' 



M=M+1 


M< MAX 



MO 

f 

COMPIITB J 

JPPCTRIM: 

PSOM AMIF(J) 


f STOP 


Plow chart for computing the MEM Spectrm 
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THE POLLOV/ING ROGTIEBS ilRE LISTED 

1. IKTERP: Interpolation by a factor L using 

M neighbouring samples 

2. SPEMM: Spectrum Estimation by Maximum Entropy 

Method (MEM) 

3* MMITTO: Direct Interpolation 
4* OPTINT: Optimimi Filter Coefficients for 
1=2 and M = 4* 

5. WGNOIZ: White Noise Series Generation - 

Uniformly Distributed or Gaussian. 



PROGRAM* 



SUBROUTINE 1 f^TERF t X f NPf H, HP f Lf Y ) 

THIS ROUTINE PERFCRNS DIGITAL INTERPOLATION BY AN INT 
X IS THE INPUT SERIES 

NP IS THE NUKBER CP INPUT DATA POINTS 

L IS THE FACTOR BY WHICH THE SAHFLING RATE IS TO INC! 
HI,J) are the IHPULSE response COEFFICIENTS WHERE I 
1.L....L ANC J RANGES FROM 1, 2,. ...HP. FOR 1=1 T! 
SAMPLES IN THE INPUT SCRIES ARE RETAINED 
MP IR the number CF input PGI^iTS TO BE USED FOR INI t 
AT PATH POINT. MP IS EVEN. FOR 1= 2» 3» •• • » MiP ^HP/2 
EITHER SIDE CF THE INTERPOLATION INSTANT ARE USED 

BINENSICN Y(L*NPWX(NPIfWlNRtNP-2l,H(L,MP) 

™ Y 1000 ) . X ( 200 ) , H ( 2 5# ) . H ( 1 0. 10 ) 


J J=MP-J + lv» 

Y(f/:)=V(MHN( 

return 

END ;L;. 






PROGRAM 


NTVLfSPECTRf MIN^ FMINs- NFREQ) 
METHOD. BURG'S ALGORITHM 
ERROR CRITERIA IS USED 
^^_JVE MODEL COEFFICIENTS AND 
each order of THE MODEL# * MAX 
rneFFlClENTS FOR THE MINIMUM F P E 
^rOEFFIClENTS ARE THEN USED IN 
AUTCCORREL at 1 ON FUtCTI ON IS TO 
F T ROUTINE FINTVL SHOULD BE SO 


SLBROUTINE SPEMEM «XtN WA,MAX»FIH * vu, t 
SPECTRUM* evaluation ’»AX1 MUM^EJTROP^ 

COUPLED: WITH AKATKEf.trTlWAJr; 

I N. THE :,FiRST'oSTAOs^TNe'‘- A utoregressive 

THE. :,FINA.I#'’‘Ptolt;tf0N:-=ER,ROR FDR 
ABE CQNP0T.E0*.'.':THE^^A -R MODEL 
ARE STORED IN AMlNIJl* THESE 
COMPUTING THE SPECTRUM, IF THE 
BE EVALUATED USING A PAD’ X 2 F 
CHOSEN THAT NFREC IS DF THE FORM 
' X IS THE INPUT SERIES ^niMTS 

NDATA is THE NUMBER OF DATA MODEL TO BE FITTED 

MAX IS THE HAXIMUM Of^DER 0 jhE SPECTRUM IS X'OMPUTEO 

FINTVL IS THE FREQ. INTERVAL AT 

TINTVL IS THE SAFPLIKG INTERVAL CF THE INPUT StKlti , 

JSy^o^TH^ Tp Uel for ;m.cH the PPE^ Him: 

FHIN IS THE HINIHUM FPE .lUrATAl SI^CTRtNFREQ) ’ " 

DIMENSION XINDATAI,B 1 (N 0 ATA 1 .B 2 ( ^ 

0 1 MENS I ON Af 12 i # A Af 12 it AMIH 1 1 i^^Tl'PIATf 

maximum ENTROPY PCWER SPECTRUM' ^ 

ALTCREGRE SSI VE MODEL '‘f' 

MINQkO-u . - rsrli-F' s'fe''.''TU; 

FMIN*0.1E + 26 ■. / 

PI = C .0 ' / ■'■'v;/- 

MAXORD=MAX • ^ ’ 

00 10 K=l»NOATA • ■■.v' 

PI*PI+XtKi=«‘XCK) 

IF(K.EQ*NCATA) GC TC 1- 
BHK) = XtKl 


aTA) 






PROGRAM 




0 

TION OF FINAL PRED CTION ERROR 


KK=IKM,NDATA 

K) 



program: 


CO 100 KQ 
AACKQ)«sA« 

CO 110: KS 


S) 

CKOSt) 




SPECTRfMI 


feSTiMATI ON 

A(Jm=AMINUM) 

FREQ=1 *0/ (TINTVL#FINTVL J 

NFREQ=INT (FREQ)+1 
DNUM=PMIN*TI NTVL 
PIE=3. 14159265 
0P0WER=2.C}*PIE*TINTVL 
AFMlN«l-0.5)/TINTVL 
00 150 MIsl, NFREC • 
ADREAL=1,C . . . 

AI31MAG«0.0 

FR€i*AFMIN+FINTVi*FtOAT|mi, 
AFRO«FRQ*DPOWER „ v.y;;;: 

DO 140 KIM«l#MlNORO' 

ANGL=AFRQ*FLOAT{K|M) 
ACREAL=ADREAL-AIKJ 
ACI*'AG=ADIMAG+ACKI 
CONTINUE 






1$ A LIBRARY SUfiRCUTlNE FCR 
’01f6N5!CN XCcO)tG(12fl2)5C{12)»Ctl2)»G£ 
CCMPLf'X r (5n5UB(A0) I Z(2aO) |TE^P0{5) 
CIFFNSICN V(80)jXX(80) 
TN-^FRPGlfiTION 
JNPsCNP+l 1/2 
jKP=''!P/2 

FVALU'vTION OF FCLRIEH TPARSFGH'F ' ; . 

DC 26 J=1?INP ■ ^ 

GG( Jt1 l^C.O ^ ^ ^ 

- DC 26 ' , . .. i. . _ . / . i,f . -im'j. .< .'Jit'' 

jk = k™j + INP ' „V, 

G(J»K)^KLaK) 1 

C Al t w A T I- V { G f I 'HP #«,f D£tf m J •; ^ ‘ 1 
CCMPUTATl-tIv CF 'i:'G'B'FE'|ip''|'£K’T| C,t 1 t.fC £ I 

DC 3 6 ' ,',t 

If t ABS'f 'Ct 

)' CCl/'PL}T^,t I CK' C'F'*;'CG£F£^^i 

'*■' 'O'C" 66' ' R *'J NP ' ' 

; ;c t : Kf K » I f R K 

‘ j' " ' ff K'F. f 4' jils F K K-'. . ■■' ,/'• j:;?/! rg'-f 
^ ' 1 6 '^p; 4 ,G |,i “ t 

■ '-'K )■■-• -M./'ifTrprsC.'c)-' 
















|'0,it6*0i(i 03 


^■l£3t»Pii:-':.oo4 . .' . 

DC life IJ*1^501 ■" ,,■ „ 

F(j j)®cppix(o*o,0e0)' .,• “ ^'3 

lll-r.-nf eX’J’FICATU J-251) /: 333331: 

DC ).?fe JK = 1,NP . *33331 

Z Z»Z77^!'FLr''.T ( JK-I KP) , , : 3 

TEMPO! Sl^C^PLXlOaOf ZZ) 

F(T J)=E (T J) + B( JK)=^CEXP(T6MPC{3) ) 
Ff!J) = TEwpc(l)/FCIJ) 

DC I4fe |J*lf500 ..,3'’*^ '',.'tf:33^ .. 

F ( T j ) :=TPMcn( 2)=?^(F !T Jj +F ( I J+1 ) ) 


» I f % * n . ^ ^ 










W 

ki * « 


:SPeCTWH"?L%trEQ«SPACED ON T« 

FRe^^Ktos’-eWEEN 


S:fil?i?®r"rE%mVR “cc^ 

'■'tilW^tiCN’ SPECTRCK) 

"‘>iii»3,«i4 

'^/FLOATiN) 

F ?=3 /’’I'fl 
S laCo4'> 

',2=G^n 

Sll=G.O 

S12 = ''*' ' Gi' 

S22*Oe‘. ill 

DC 6nf,0 T=1,N 11 

TfcH = criS(Fl*CLCAT ? I ) ) ’^SPECTP 1 1 V..,.,^-';; 

•:]:rC] + TFF I . 

T 6M»crjst 

^ ^ ^ ^ , . wi*is'L citif 1 -I r mi ..-;i: i ; - 



■wm- 




DO 70 

DO 60 I *1 tWZ ■ ' : 

^"X2«xtKiN2)-':=cexp'( 






-- vj='y-y ^^ -■'.; 


i-lklliWf.tOA'TIH) ip ’’•' 

H ‘ r '''^. / «■''. kvik '‘■^'’?1, . ‘f 




',^<5rs. i» •‘^ 


k=k+h: 


''f|N»N 
f)0 10 l»lfH 
NN=NN/2 
K=M+1 

IFfNN„LT.2) GO TO 20 
COMTJ'IUE 

I F( IMVFFT«NE.I) go to 40 
DC 30 1=1 »N 
X(I)=CONJG(X(l)) 

N2=M/? 



FRCGRAf-' 


GC TC eO 

'XIK)=*Xtn^’ 

XU )*X 3 
‘'r'XCNTTMtE 

S'i!'''.,; I f(TNVFFT.KE«1 5 GC TC. ICC 

' DC ‘SC K = 1^K 

9C x <K)=f.CNJG(X(K) ) 

ICC RETLRM 

END 

FIA'CTT CN T SCR'^'l ( s! 

J2*J 

ISCRWL=C 
DC nn 1 = 1, V 

' J 1 ® J 2 / 2 

lie -■ 42’Jl • AC 


£ND 




'• ' '' -■ '' 






.; ' I,.'*' ‘ t', i ' * •* ' i\'‘ T -*4' 


RTJHGP.AS^ 




l'&WO /'B w "c 
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